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We show that, along the critical isobar, the complete scaling results in a unique leading scaling 
qualitatively distinct to that arising from the simple and the revised scalings. This is verified by 
a complete-field finite-time scaling theory, which combines the complete scaling with finite-time 
scaling, and its application to the molecular dynamics simulations of the vapor-liquid critical point 
of a three-dimensional one-component Lennard-Jones fluid in an isobaric-isothermal ensemble with 
linear heating or cooling. Both the static and the dynamic critical exponents as well as the critical 
parameters can be estimated without a priori knowledge of the universality class. The results agree 
with extant values and thus show the necessity of the complete scaling to the leading asymptotic 
behavior along the critical isobar even for the LJ fluid whose asymmetry is thought to be weak. 
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Fluid criticality still attracts great attention though it 
was first recorded nearly one and a half century ago [1]. 

A primary concern is the lack of symmetry between the 
vapor and the liquid phases for usual fluids and the rel¬ 
evant complicated field mixing [2]. 

In the “simple scaling”, no mixing is needed and the 
reduced temperature r = (T — T^jT^ and the chemi¬ 
cal potential /t = (/i — Tc)I^p,Tc correspond directly to 
the thermal and the ordering fields in the Ising model, 
where fee is the Boltzmann constant and the subscript 
c stands for the values at the critical point hereafter. 
This leads to the usual singularity of the order param¬ 
eter A/5 = {p — Pc)/Pc ~ \t\^ with its critical exponent 
/3, where p is the density. Yet, along the vapor-liquid 
coexistence curve this leading behavior cancels and the 
resultant diameter pd = (pvap + Piiq)/2 = Pc is sym¬ 
metric, different from the asymmetry in the empirical 
law of rectilinear diameter pd oc |r| [3]. In the “revised 
scaling” [4], mixing of the two physical fields gives rise 
to a subleading term in Ap, a term which accounts for 
the singularity in pa [5, 6], where a is the crit¬ 

ical exponent of the isochoric heat capacity. Recently, 
a “complete scaling” has been proposed [7-10] in which 
the pressure P is further mixed into the scaling fields in 
order to account for the Yang-Yang (YY) anomaly [II]. 
This results in yet another term in A/3 and hence an 
additional [rp^ singularity in pd, which is dominant as 
2/3 < 1 — a normally. The complete scaling is supported 
by data on vapor-liquid coexistence and heat capacity in 
highly asymmetric fluids [12], and is found to be equiv¬ 
alent to a field-theoretic treatment of asymmetric fluid 
criticality up to independent fifth-order terms [13]. Some 
experiments from binary fluids also support the 2/9 sin¬ 
gularity [14-18] and the YY anomaly [19, 20]. However, 
experimental evidences from specific heat [21] for the YY 
anomaly may stem from small traces of impurities [22], 
and the two singular terms may compensate in some flu¬ 


ids, resulting in an almost rectilinear diameter [12]. As 
the leading singularity of A/5 remains intact, a question 
then arises as to whether it could be changed so that the 
complete scaling could be more detectable. 

Another concern is computer simulations. Although 
simulational estimates of the critical properties of lattice 
systems such as Ising models and the like has reached 
a high level of accuracy, that of off-lattice systems such 
as fluids is less satisfactory. Besides the asymmetry, an¬ 
other main difficulty arises from the divergent correlation 
length ^ for fluctuations at critical points. Usual methods 
such as direct interfacial [23, 24], Gibbs ensemble Monte 
Carlo (MC) [25], isobaric-isothermal ensembles [NPT, 
N is particle number) [26, 27] or canonical {NVT, V is 
volume) [28, 29] ensembles plus test particle [30] usually 
work away from the critical point due to strong fluctu¬ 
ations there and the accompanying extrapolation to de¬ 
termine critical properties is questionable [31, 32]. 

Near to critical points, finite-size scaling (FSS) [33- 
35] is crucial due to the divergent MC simulations in 
grand canonical [pVT) ensembles, amenable to an FSS 
analysis, are thus usually employed [31]. By contrast, 
NPT ensembles need a revised FSS and are deemed to be 
computationally less efficient [36]. A mixed-field FSS the¬ 
ory [31, 37, 38] based on the revised scaling for estimat¬ 
ing critical properites of fluids has been developed using 
pVT ensemble MC simulations with histogram reweight¬ 
ing [39]. Yet, this method is “biased” in that a pri¬ 
ori knowledge of the Ising universality class ought to be 
given [40] . An unbiased recursive algorithm based on the 
complete scaling using the FSS of cumulants [41] has also 
been developed again with MC simulations in pVT en¬ 
sembles [42-44]. However, the method demands a signif¬ 
icant amount of high quality precise data and adjusting 
parameters to be determined to exact or approximated 
scaling functions. In addition, dynamic critical behavior 
has yet to be incorporated. 
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When dynamics is considered one encounters the diver¬ 
gent relaxation time teq that brings about critical slowing 
down as critical points are approached. Also, critical dy¬ 
namics of fluids is complicated by hydrodynamic modes 
and is described by model H [45] . Consequently, although 
there exist consistent theoretical [46-48] and experimen¬ 
tal [49, 50] results, molecular dynamics (MD) simula¬ 
tions for the critical dynamics of a Lennard-Jones (LJ) 
fluid [51] and a Widom-Rowlinson mixture [52] appeared 
merely about a decade ago. The former used a specially 
designed, near equilibrium initial condition in an NVT 
ensemble but knowledge of equilibrium critical properties 
is a prerequisite. The latter employed an event-driven 
MD algorithm to expedite the dynamics [53]. However, 
the value of the dynamic critical exponent z obtained 
disagrees with the accepted ones [54]. A correction by 
background remedies the discrepancy but theoretical val¬ 
ues have to be invoked using microcanonical ensemble 
MD simulations [55]. In the latter two pieces of work, 
semi-fiVT ensembles [56] are employed to provide both 
equilibrium properties and initial conditions for the sub¬ 
sequent MD calculations. One sees that MD simulations 
alone could not yet produce all critical properties though 
they simulate the real dynamic process. 

Here, we show that the complete scaling can result in 
a qualitatively distinct leading behavior to the other two 
scalings along the critical isobar for Ap itself and other 
thermodynamic functions besides pd- As P is mixed into 
the scaling fields in the complete scaling, if we fix it and 
approach the critical point along the critical isobar, both 
scaling fields are affected, in sharp contrast to the other 
two scalings. To show this, we explore the time aspect 
and develop a theory of “complete-field” finite-time scal¬ 
ing by combining finite-time scaling (FTS) [57, 58] with 
the complete scaling. Moreover, we are able to effect the 
theory with MD simulations of a three-dimensional pure 
LJ fluid. NPT ensembles are employed though T is var¬ 
ied linearly with time to realize FTS, a situation which 
can be realized in experiments. Our numerical results 
support to the leading asymptotic behavior the complete 
scaling, which is thus needed along the unique critical 
isobar even for the LJ fluid whose asymmetry is believed 
to be weak. We can also determine both the static and 
dynamic critical exponents as well as the critical param¬ 
eters with the dynamic NPT ensembles alone. 

FTS has been applied effectively both to classical [57- 
62] and quantum phase transitions in lattice models [63- 
65] . Analogue to FSS that circumvents the problem that 
^ is longer than the lattice size L, FTS overcomes the 
critical slowing down of a divergent teq by devising a con¬ 
trollable time scale tdr- This tdr is obtained by linearly 
varying T at a constant sweep rate R to drive a system 
through its critical point. Although tdr is shorter than teq 
and thus the system falls out of equilibrium in the FTS 
regime, the FTS enables one to extract both static and 
dynamic critical properties similar to its spatial counter¬ 


part. Moreover, in the FTS regime, the effective length 
scale arising from the driving is shorter than L and thus 
finite-size effects are only subsidiary and negligible. 

To see this, note that the time evolution of the char¬ 
acteristic thermodynamic function in an NPT ensemble 
under a coarse-graining of a factor b is 

p = ,th-^), ( 1 ) 

where p, f, and p are the scaling fields corresponding 
to p, T, and p = [P — Pc)/P c^bTc, respectively, and <5 
and v are the critical exponents of the critical isotherm 
and respectively. For simplicity of presentation, we 
shall use identical symbols for all the three scalings and 
differentiate the latter by their contents. Thus, in the 
simple scaling, all the three scaling fields are just the 
physical fields themselves, while in the complete scaling 
in the linear approximation [7, 8, 12], 

p = p-koT-lop, f = T-lip-jip, p = p-kiT-j2p, 

( 2 ) 

where all coefficients are constants. In the revised scaling, 
ko, Iq, ji, and j 2 all equal zero. In order to avoid reaching 
an apparently different universality class, we have asso¬ 
ciated the field variable p and the function p in the NPT 
ensemble with the exponents 2 —a and /JJ, respectively, a 
reversal of the usual assignment in magnetic systems [66] . 

In the simple scaling, Eq. (1) describes the usual crit¬ 
ical dynamics; but it is valid even if t = Rt, viz., time 
dependent, which has been established by a dynamic 
renormalization-group (RG) theory [57, 58, 67]. Replac¬ 
ing T with R and choosing b such that RN'^ is a constant, 
one arrives at an FTS form 

where fi (i an integer) is a scaling function and rx is 
the RG eigenvalue of R. We have neglected dimensional 
factors and the possible difference in fi between the two 
phases hereafter for simplicity, tt = zj-l/zz from t = Rt 
and its coarse-grained version [67]. Equation (3) shows 
that the driving imposes a timescale tdr on the 

evolution of the system. In the FTS regime, ^ 1, 

or, back to r, teq ^ ^dr as expected. If one 

includes L in Eq. (1), one will see that the ETS regime is 
given additionally hy L ^ R~^At ^ which is an effective 
driven length scale. So, a series of driving circumvents 
both the divergent length and the divergent time. 

Next we turn to the critical isobar and develop the 
complete-field FTS theory by incorporating the complete 
scaling in FTS. Note that all scaling fields vary with r 
and the method that leads to Eq. (3) does not work. So, 
we need to find the leading behaviors of all the scaling 
fields with r at P = Pc- This can be done by expanding 
the equilibrium scaling form, 

P=\f\^^h{P/\T\^-n. 


( 4 ) 
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which results from Eq. (1), in p for vanishingly small f 
and using ( 2 ), resulting to the leading order in [ 8 ] 

f ^ r, A ~ ~ |f |(2-a)//3<5_ 

Since Ap = {dp/dfi)r — 1, it can be related to the scaled 
density {dp/dfi)f through Eq. (2). So, 

Ap = Ap{n^l'' ( 6 ) 

similar to Eq. (1) to the leading asymptotic term [7, 8 , 
12]. Substituting Eq. (5) into ( 6 ), choosing a scale such 
that tb~^ is a constant, and expressing t in terms of r 
and i?, we finally obtain the leading singularity 

Ap(T,i?) = 

_ (7) 

where r = z + /3S/v and in comparison to the neglected 
subleading terms we have neglected in the second line 
a correction to scaling since fiS/rv — I/ttv > 0. Note 
that r is the rate exponent associated with varying the 
ordering field in magnetic systems [57, 58, 67]. This can 
also be seen if we consider instead a p,VT ensemble and 
vary p, in the simple scaling (see below). 

The FTS form (7) is qualitatively distinct from those 
for the simple and the revised scalings in the same pro¬ 
cess. In the simple scaling, Eq. (3) results obviously in an 
FTS along the critical isobar characterized by the ther¬ 
mal rate exponent rr instead of the field rate exponent r. 
In the revised scaling, p = p = 0 then leads to p ^ Jf]^^ 
from Eq. (4). This then gives rise to p, ^ t and so r ^ r 
from Eq. (2). So, from Eq. ( 6 ), varying r can only result 
in an FTS form similar to the simple scaling. 

We now show that the critical isobar with its result 
is unique in the sense that along other paths such as 
the coexistence curve or the critical isokyme, the leading 
behavior of the complete scaling cannot be uniquely dis¬ 
tinguished from the others. Along the coexistence curve, 
chosen as p = 0 [12], p ~ from Eq. (4). As a re¬ 
sult, in Eq. ( 6 ), is a function of fb^^‘' and thus 

only tt matters. So do the other two scalings. Along 
the critical isokyme, viz. fixing p to pc and varying T, 
Eq. (5) is also valid [ 8 ] and thus Eq. (7) applies as well. 
This is also true for the revised scaling, though in the 
simple scaling, only rx appears as both paths coincide. 

Note that the above conclusions do not depend on the 
ensemble used. In a pVT ensemble, p 6 ( 2 -a)/i/ ^ 0 ^ 

is replaced by pb^^^^. Along the critical isobar, for the 
complete scaling, the first two expressions in (5), valid in 
this case too, leads again to (7). For the revised scaling, 
again p ~ \t\^^ and f ^ t. Accordingly, both pb^^^'^ 
and fb^/'^ are reduced to rb^/'^ and thus the FTS form 
involves rx instead of r, similar to the simple scaling. 
Note also that the results are not due to FTS. In equi¬ 
librium, one can see that the complete scaling leads to 
a leading behavior of Ap ^ [ 8 ], whereas the other 


two result in Jr]^, reflecting both the mixing of P and 
the dominating of the field. 

Having known the FTS form of Ap, those for other 
observables can be similarly derived. For example, the 
isobaric specific heat Cp and the reduced entropy density 
AS = {S — Sf)/Pcks along the critical isobar, the for¬ 
mer calculated from the fluctuations of enthalpy and the 
latter from the integration of the former, behave as 

Cp(T, R) - Cpo 2=^ (8) 

AS{t, R) ~ (9) 

respectively, to the leading singularity, where Cpo is a 
background contribution. 

From the scaling forms, one can then determine the 
critical properties using the FTS methods [57, 58]. We 
utilize p and S instead of the peaks in Cp as they are 
smoother. If Pc has been known, one can use Eqs. (7) and 
(9) and a method of iterations to estimate the parameters 
involved. In particular, given an initial Tc, one obtains 
p at different sweep rates at this Tc- Then pc and fjrv 
can be extracted from 

p = PcParR^I’^\ ( 10 ) 

according to Eq. (7), where {i an integer) is a constant. 
At this pci the temperature T at each rate can also be 
identified to find a new Tc and fb/rv from 

r = Tc-b (11) 

again from Eq. (7). Several iterations suffice to reach 
convergent Tc and pc and thus j35 /rv and f jrv. A similar 
procedure can be conducted for S and thus Sc^ {l—a)/rv, 
and again Tc and fbjrv can be obtained, the latter two 
parameters can serve as a consistency check. 

We have yet to determine Pc- Note that at a criti¬ 
cal point, Tc obtained from the peak temperatures of Cp 
using Eq. (11) during heating and cooling is identical; 
however, it is different below Tc [68, 69]. This can be 
employed to estimate Pc [60]. 

Our model system contains atoms interacting via the 
LJ potential fix) = 4e[(o-/x)^^ — (cr/x)®], where x is the 
distance between two atoms and e and a specify the units 
of energy and length, respectively. Thus, we use reduced 
quantities denoted by asterisks such as energy E* = E/e, 
time t* = t/(m(7^/e)^/^, L* = L/a, T* = ksT/e, and 
P* = Pa^/e, where m is the mass of an atom. The full 
LJ potential is used with a cutoff at 3.0cr. Berendsen’s 
method is applied to couple T* and P* of the system to 
the bath [70]. We employ the leap-frog algorithm [71] 
with a time step 5t* = 0.001. The system contains 150 
atoms. At each given P*, the system is heated or cooled 
from an equilibrium state to the other state by linearly 
varying T* at a series of prescribed R* from 5 x 10“^° to 
1 X 10“®. Larger rates result in undistinguishable evolu¬ 
tion but lower rates are time consuming and may cross 
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FIG. 1. (Color online) T* extracted from Eq. (11) at various 
P* from heating and cooling. Lines are a guide to the eye. 

over to an FSS regime [58, 62]. Several thousand sam¬ 
ples with different initial states are run for average. At 
P* = 0.140, only heating is conducted and the sample 
size is 10 , 000 . 

Figure 1 shows the estimation of P*. The difference 
between T* from heating and cooling diminishes as P* 
increases and vanishes within the errorbars between P* = 
0.14 and 0.145. Note that above P* no transition exists 
at all and the estimated T* is at best a crossover. We 
thus estimate P* = 0.142(3) with T* = 1.349(2). 

Having determined P*, we then apply the method 
of iterations to find other parameters. The results are 
pI = 0.308(1), P/rv = 0.094(1), T* = 1.348(1), and 
f35/rv = 0.448(5) from p* and S* = 0.0076(1), (1 — 
a)/ri^ = 0.253(1), T* = 1.351(2), and (3S/ri^ = 0.445(2) 
from S*. Consistencies in T* and jSS/rv are manifest. 
These results yield reasonably data collapses according 
to Eqs. (7) and (9) as seen in Figs. 2 and 3 for the aver¬ 
aged density and entropy, respectively, noting that only 
the leading contributions have been considered. Using 
7 = /3((5 — 1), we can also plot data collapses for C* ac¬ 
cording to Eq. ( 8 ) for a further check. The background 
is set at C*q = 0.083, the common C* of different R at 
the lowest T*, as C*q must change little in the vicinity 
of T*. The rescaled curves overlap quite well as shown 
in Fig. 4 and corroborate our results. 

All individual critical exponents can be derived from 
the three known ratios using the scaling laws a+l3+/3S = 
2 and a = 2 — 3;/. The results are <5 = 4.74(7), 7 = 




FIG. 3. (Color online) Reduced entropy density versus tem¬ 
perature (a) and its rescaled (b) for various sweep rates at 
P*. Inset: The entire rescaled curves. 


1.226(28), /3 = 0.328(7), = 0.627(5), a = 0.119(16), all 

in good agreement with 4.789(2), 1.233(22), 0.3265(3), 
0.6301(4), and 0.110(1), respectively, of the Ising univer¬ 
sality class [72], and z = 3.078(40), also consistent with 
theoretical [46-48] , experimental [49, 50] , and the recent 
simulational [51] results. In addition, r = 5.556(52) and 
rx = 4.673(42), quite different. This indicates that ab¬ 
sence of the complete scaling must lead to far poorer 
results if only the leading behavior is considered. 

The nonuniversal critical parameters obtained agree 
with p* = 0.304(6) [73], T* = 1.35 [73], and P* = 
0.147(2) [74] individually albeit differing somehow from 
pI = 0.316(2) [27, 73-78], P* = 0.13(1) and T* = 1.32(1) 
in the thermodynamic limit [74-79]. Several factors are 
of concern. One is the precision of Pc- Within the theory, 
one can show, however, that in the FTS regime, small de¬ 
viations of Pc do not modify the results, while too close 
to Tc due to too small R will lead to a crossover to a 
regime entirely determined by the deviations, a regime 
in which the i?-dependence is only a correction. This is 
similar to the case where there is a residue magnetic field 
in magnetic systems [80]. Others include contributions 
from the subleading terms and the small size [58, 62]. 
Further studies are needed to clarify their effects. How¬ 
ever, they should not change our main conclusion as the 
relative exponents are divided by the large rate exponent 
in FTS and are thus small [60, 61]. 



FIG. 2. (Golor online) Reduced density versus temperature FIG. 4. (Color online) Reduced isobaric specific heat versus 
(a) and its rescaled (b) for various sweep rates listed at P*. temperature (a) and its rescaled (b) for various sweep rates 
Inset: The entire rescaled curves. at P*. Inset: The entire rescaled curves. 
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Summarizing, we have shown that along the critical 
isobar, the complete scaling results in a unique leadiirg 
scaling different from the simple and revised scaling. We 
have also developed a complete-field FTS theory for fluid 
criticality and a corresponding method to estimate both 
the static and the dynamic critical exponents as well as 
critical parameters using molecular dynamics simulations 
in an isobaric-isothermal ensemble alone. 

We greatly appreciate J. V. Sengers, M. A. Anisimov, 
A. Angell, and S. K. Das for their information and dis¬ 
cussions. This work was supported by the NNSF of PRC 
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